﻿extern alias alglibnet2;

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;


namespace ClusterAggregation.Clusterers
{
    using Datum;
    using DataSets;
    /**
     * Spectral Clusterer implementation
     * has 2 modes.
     * mode 1 is without laplacian (smoothing parameter = 0)
     * the other one is with laplacian.
     * Take note that the smoothing parameter in code is actually 1/smooth_param for better accuracy, and as such it should belong to (0,1].
     *
     */
    public class CSpectralClusterer:IClusterer
    {
        /** k-means needed to finalize the clustering */
        private CKMeansImpl m_kmeansClusterer;
        /** default distance function (for the case we receive null) */
        private ISimilarity m_distFunc;
        /** smoothing parameter */
        private double m_smoothParam;
        /** this param, when using laplacian, decides if we will normalize the similarity matrix */
        private bool m_normalize;
        /** a default clusterer, in case AlgLib doesn't solve eigenvectors */
        private static IClusterer m_sDefaultClusterer = new CSpectralClusterer(null,0,false);
        /** a secondary default clusterer, in case AlgLib doesn't solve eigenvectors */
        private static IClusterer m_sSecondaryDefaultClusterer = new CPamImpl();
        /**
         * basic constructor.
         * @param distFunc (ISimilarity) [IN] default distance function.
         * @param smooth (double) [IN] sets the smoothing param to be -2 * smooth ^ 2 
         * @param normalize (bool) [IN] this param, when using laplacian, decides if we will normalize the similarity matri
         */
        public CSpectralClusterer(ISimilarity distFunc = null, double smooth = 0.7, bool normalize = true)
        {
            m_kmeansClusterer = new CKMeansImpl();
            m_distFunc = distFunc;

            if (distFunc == null)
                m_distFunc = new CEuclideanDistanceSimilarityFunction();

            m_smoothParam =  ( smooth * smooth * (- 2 ));

            m_normalize = normalize;
        }

        public CPartition cluster(AData[] arg, ISimilarity distFunc, int k)
        {
            lock(this)
            {
                DateTime dt = DateTime.Now;
                double[,] paramMat = distFunc.getFullMatrix(arg);
                double[] wr = new double[arg.Length]; // real part of eigenvalues  (no imaginary part since we send a symetric matrix
                double[] dataRank = new double[arg.Length];
                double[,] eigenVectors = new double[arg.Length, arg.Length];
                int breaker;
                double maxjump = 0;
                int ind;

                // get ranks
                for (int i = 0; i < arg.Length; i++)
                {
                    dataRank[i] = 0;
                    for (int j = 0; j < arg.Length; j++)
                    {
                        dataRank[i] += paramMat[i, j];
                    }
                }

             
                //init param mat
                if (m_smoothParam == 0)
                    m_smoothParam = 0 ; // nothing happens at all

                else
                {
                    //calculate diagonal values of D ^ -0.5
                    for (int i = 0; i < arg.Length; i++)
                    {
                        dataRank[i] = Math.Pow(dataRank[i], -0.5);
                    }

                    //// adjust rank matrix a little bit;

                    if (m_normalize) // normalize paramMat
                    {
                        // begin by calculating sum of each row
                        double[] sum = new double[arg.Length];
                        for (int i = 0; i < arg.Length; i++)
                        {
                            sum[i] = 0;

                            for (int j = 0; j < arg.Length; j++)
                            {
                                sum[i] += paramMat[i, j];
                            }
                        } // end of calculate sum of rows
                        for (int i = 0; i < arg.Length; i++)
                        {
                            for (int j = 0; j < arg.Length; j++)
                            {
                                paramMat[i, j] /= sum[i > j ? i : j];
                            }
                        }
                        // end of normalize paramMat 
                    }

                    // calculate W
                    for (int i = 0; i < arg.Length; i++)
                    {
                        for (int j = 0; j < arg.Length; j++)
                        {
                            paramMat[i, j] = Math.Pow(Math.E, paramMat[i, j] * m_smoothParam);
                        }
                        //paramMat[i, i] = 0; - similarity between i and i is 1! don't override
                    }


 
                    // calculate W*D
                    for (int i = 0; i < arg.Length; i++)
                    {
                        for (int j = 0; j < arg.Length; j++)
                        {
                            paramMat[i, j] = paramMat[i,j] * dataRank[j]; 
                            // Line I, row J, in row J only the J'th number is not 0
                        }
                    }
                    // calcuulate D*(W*D)
                    for (int i = 0; i < arg.Length; i++)
                    {
                        for (int j = 0; j < arg.Length; j++)
                        {
                            paramMat[i, j] = paramMat[i, j] * dataRank[i]; 
                            // more understandable if you're doing it on a piece of paper
                        }
                    }// 
                }

                // calc eigen vectors and values              
                bool resp = alglibnet2.alglib.evd.smatrixevd (paramMat, arg.Length , 1/* eigenvectors needed */,true,ref wr,ref eigenVectors);

                breaker=ind = arg.Length - 2;
                
                while (ind >0 && wr[ind] > 0)
                {
                    if (maxjump < Math.Abs(wr[ind] - wr[ind + 1]))
                    {
                        maxjump = Math.Abs(wr[ind] - wr[ind + 1]);
                        breaker = ind;
                    }
                    ind--;
                }
                int top = arg.Length - 1;

                AData[] toCluster = new AData[arg.Length];

                breaker = top - breaker + 1; // less calculations later

                for (int i = 0; i < arg.Length; i++)
                {
                    toCluster[i] = new AData(breaker); // allocate coordinates

                    toCluster[i].id = i;
                    for (int j = 0; j < breaker; j++)
                        if (double.IsNaN(eigenVectors[i,j]))
                        {
                            if (m_smoothParam == 0)
                                return m_sSecondaryDefaultClusterer.cluster(arg, distFunc, k);
                            else
                                return m_sDefaultClusterer.cluster(arg, distFunc, k);
                        }
                        else
                            toCluster[i].coordinates[j] = eigenVectors[i,j];  
                }

                CPartition tmpPartition = m_kmeansClusterer.cluster(toCluster,m_distFunc,k);
                CPartition finalPartition = new CPartition();
                CCluster tmpCluster;
                foreach (CCluster clstr in tmpPartition._clusters)
                {
                    tmpCluster = new CCluster();
                    for (int i = 0; i < clstr._data.Count; i++)
                    {
                        // clstr._data[i].id  is the id of the data we want to add
                        // we add from arg
                        tmpCluster._data.Add(arg[clstr._data[i].id]); 
                    }
                    finalPartition._clusters.Add(tmpCluster);
                } // end of foreach
                DateTime dt2 = DateTime.Now;
                int hour = dt2.Hour - dt.Hour;
                int min = dt2.Minute - dt.Minute;
                int sec = dt2.Second - dt.Second;
                finalPartition.name = "Spectral Clustering";
                return finalPartition;
            } // end of lock
        }//end of func

    }
}
